/* Poisson equation solver using red-black Gauss-Seidel successive over-relaxation iteration. */

# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

#ifdef _OMP_SOLVE_POISSON_
# include <omp.h>
#endif

/*
# include <time.h>
# include <sys/timeb.h>
*/

/* This structure is used to detect most sensitive cell in an itteration. */
struct dpmax
    {
     short im, jm, km;
     double dpm;
    };

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK);

extern long int solve_poisson(struct grid *g, const double ***d, unsigned short iRK)
  {
   int cnt;
   short i, j, k, im, jm, km, stop_cnt, tid, nTRD;
   double dpm, dp, ***p, *r, *th, *z, prf;
   char flag;
   struct dpmax *tdp;

// struct timeb t1,t2;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'solve_poisson()'...\n");
#  endif

   r = g->r; z = g->z; th = g->th;

   prf = PRF_INIT; stop_cnt = 0;

   nTRD = 1;
#  ifdef _OMP_SOLVE_POISSON_
   nTRD = g->nTRD;
#  endif

   tdp = (struct dpmax *)malloc(sizeof(struct dpmax)*nTRD);

   p = g->p;
   for(cnt = 1; cnt <= POISS_CNT_MAX; cnt++)
      {
#      ifdef _OMP_SOLVE_POISSON_
       #pragma omp parallel private(i,j,k,dp,tid)
#      endif
          {
        /* Initializing elements of the structure dpmax. */
#          ifdef _OMP_SOLVE_POISSON_
           tid = omp_get_thread_num();
#          else
           tid = 0;
#          endif
           tdp[tid].dpm = (double)0;
           tdp[tid].im = -1;
           tdp[tid].jm = -1;
           tdp[tid].km = -1;

#          ifdef _OMP_SOLVE_POISSON_
           #pragma omp for schedule(guided)
#          endif
           for(k = NzG; k < NzG+g->nz; k++) for(i = NrG; i < NrG+g->nr; i++) for(j = NthG+(k+i)%2; j < NthG+g->nth; j += 2)
              {
               dp = P_ijk_NEW - p[i][j][k];

               if((dp<0 ? -dp:dp) > (tdp[tid].dpm<0 ? -tdp[tid].dpm:tdp[tid].dpm))
                  {
                   tdp[tid].dpm = dp;
                   tdp[tid].im = i;
                   tdp[tid].jm = j;
                   tdp[tid].km = k;
                  }
               p[i][j][k] += prf*dp;
              }

#          ifdef _OMP_SOLVE_POISSON_
           #pragma omp for schedule(guided)
#          endif
           for(k = NzG; k < NzG+g->nz; k++) for(i = NrG; i < NrG+g->nr; i++) for(j = NthG+(k+i+1)%2; j < NthG+g->nth; j += 2)
              {
               dp = P_ijk_NEW - p[i][j][k];

               if((dp<0 ? -dp:dp) > (tdp[tid].dpm<0 ? -tdp[tid].dpm:tdp[tid].dpm))
                  {
                   tdp[tid].dpm = dp;
                   tdp[tid].im = i;
                   tdp[tid].jm = j;
                   tdp[tid].km = k;
                  }
               p[i][j][k] += prf*dp;
              }
          }

    /* Calculating the absolute maximum. */
       dpm = tdp[0].dpm;
       i = 0;

       for(tid = 1; tid < nTRD; tid++)
         if( (tdp[tid].dpm < 0 ? -tdp[tid].dpm:tdp[tid].dpm) > (dpm < 0 ? -dpm:dpm) )
             {
              dpm = tdp[tid].dpm;
              i = tid;
             }
       im = tdp[i].im;
       jm = tdp[i].jm;
       km = tdp[i].km;

       press_bc(g, 0, iRK);

       if( (cnt == 1) || (cnt%DISP_FREQ == 0) )
            {
             fprintf(g->dptr, "=> %d: dpm = %E(%d, %d, %d)\n", cnt, dpm, im, jm, km);
             fflush(g->dptr);
            }

       if( ((dpm < 0) ? -dpm:dpm) <= PRESS_DELTA) ++stop_cnt;
       else stop_cnt = 0;

       if(stop_cnt == POISS_CNT_MIN) break;
      }

   free(tdp);

   flag = (stop_cnt == POISS_CNT_MIN) ? 'Y':'N';

   fprintf(g->dptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c\n", cnt, dpm, im, jm, km, flag);
   fflush(g->dptr);

   fprintf(g->tptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c; ", cnt, dpm, im, jm, km, flag);
   fflush(g->tptr);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'solve_poisson()' ends.\n");
#  endif

   return cnt;
  }
